function [wstar2,cumcomp,cumulsum,meanacceptedw,expectedw,proboffer,probs2,rej_atleast1,rej,fdist,fval,wgriduse,ecost2,expectedoff,cstar2,jfr] = ...
    modelsolve_endosearch(iota,b,zeta,mu0,gam,mustar,sigstar,scost,lam,bet,capT,capT2,Nw,Nwlong,cumcomp_dat,beliefs_dat,proboffer_dat,cumsum_dat)

warning('off')
options = optimset('MaxFunEvals',1e3,'TolX',1e-7,'MaxIter',1e3,'TolFun',1e-8) ;

% beliefs over time %
mu(:,1)           = mu0.*(exp(-gam.*(linspace(0,capT2-1,capT2))'))+mustar.*(1-(exp(-gam.*(linspace(0,capT2-1,capT2))'))) ;


% wage grid %
wlb            = .001 ;
wub            = logninv(.99999,max(max(mu)),sigstar(1)) ;
[wgrid(:,1),~] = chebpts(Nw,[wlb wub]);
[wgriduse(:,1), wgridwgts(:,1)] = chebpts(Nwlong,[wlb wub]) ;  


wub2            = logninv(.99999,mustar,sigstar(1)) ;
[wgrid2(:,1),~] = chebpts(Nw,[wlb wub]);
[wgriduse2(:,1), wgridwgts2(:,1)] = chebpts(Nwlong,[wlb wub2]) ;  


% initial guesses %
W                           = zeros(Nw,capT2) ; % wages by job type time (not a function of beliefs)
U                           = zeros(capT,capT2) ; % time by beliefs 
U(capT2,1:capT2)            = (((b.*(1-zeta)).^(1-iota)-1)./(1-iota)) ;
W(:,capT2)                  = (((wgrid.^(1-iota)-1)./(1-iota)))./(1-bet) ;
wstar(capT2,1:capT2)        = NaN ;
probs(capT2,1:capT2)        = NaN ;

for m = 1:capT2 % for each belief 
    for t = capT2:-1:1
        if t==capT2
                tol = 100 ;
                while tol>1e-7
                    [w1]                            = ndgrid(wgrid(:,1)) ;
                    Wuse                            = griddedInterpolant(w1,W(:,t),'spline','spline') ;
                    [hi1,~]                         = find(min(((W(:,t)-U(t,m)).^2))==(W(:,t)-U(t,m)).^2,1,'first') ;
                    try
                        [ystar] = fminsearch(@(y) bisect_wage(y,Wuse,U(t,m)),log(b),options) ;
                    catch
                        display('problem')
                    end
                    wstar(t,m)                      = exp(ystar) ;
                    int                             = bet.*integral(@(w) lognpdf(w,mu(m,1),sigstar).*(Wuse(w)-U(t,m)),wstar(t,m),Inf);
                    cstar(t,m)                      = int*lam ;
                    probs(t,m)                      = 1-exp(-scost*(lam*(int))) ;
                    ecost(t,m)                      = expected_cost(lam*int,scost) ;
                    U1    = probs(t,m)*(-ecost(t,m)+((b*(1-zeta))^(1-iota)-1)./(1-iota) ...
                                    +lam.*bet.*integral(@(w) lognpdf(w,mu(m,1),sigstar).*Wuse(w),wstar(t,m),Inf)+lam.*bet.*U(t,m).*logncdf(wstar(t,m),mu(m,1),sigstar)+ ...
                                    bet.*(1-lam).*U(t,m))+...
                                    (1-probs(t,m))*(((b*(1-zeta))^(1-iota)-1)./(1-iota)+bet.*U(t,m));

                    tol = abs(U1-U(t,m)) ;
                    U(t,m) = U1 ;
                end 
        elseif t>=capT
            U(t,m) = U(capT2,m) ;
            wstar(t,m) = wstar(capT2,m) ;
            probs(t,m) = probs(capT2,m) ;
            ecost(t,m) = ecost(capT2,m) ;
            cstar(t,m) = cstar(capT2,m) ;
            W(:,t)     = W(:,capT2) ;
        else
            W(:,t)         =  ((b).^(1-iota)-1)./(1-iota)+bet.*W(:,t+1) ;
            [w1]           = ndgrid(wgrid(:,1)) ;
            Wuse           = griddedInterpolant(w1,W(:,t+1),'spline','spline') ;
            [hi1,~]        = find(min(((W(:,t+1)-U(t+1,m)).^2))==(W(:,t+1)-U(t+1,m)).^2,1,'first') ;

            [ystar]     = fminsearch(@(y) bisect_wage(y,Wuse,U(t+1,m)),log(b),options) ;
            wstar(t,m)     = exp(ystar) ;
            int            = bet.*integral(@(w) lognpdf(w,mu(m,1),sigstar).*(Wuse(w)-U(t+1,m)),wstar(t,m),Inf) ;
            cstar(t,m)                      = int*lam ;
            probs(t,m)                      = 1-exp(-scost*(lam*(int))) ;
            ecost(t,m)                      = expected_cost(lam*int,scost) ;

            U(t,m)  = probs(t,m)*(-ecost(t,m)+((b*(1-zeta))^(1-iota)-1)./(1-iota) ...
                                    +lam.*bet.*integral(@(w) lognpdf(w,mu(m,1),sigstar).*Wuse(w),wstar(t,m),Inf)+lam.*bet.*U(t+1,m).*logncdf(wstar(t,m),mu(m,1),sigstar)+ ...
                                    bet.*(1-lam).*U(t+1,m))+...
                                    (1-probs(t,m))*(((b*(1-zeta))^(1-iota)-1)./(1-iota)+bet.*U(t+1,m));     
        end
               
    end
end

wstar2 = diag(wstar) ;
probs2 = diag(probs) ;
ecost2 = diag(ecost) ;
cstar2 = diag(cstar) ;
wint = @(x) x.*lognpdf(x,mustar,sigstar);
for t = 1:capT2
    rej(t,1)           = probs(t,t)*lam*(logncdf(wstar(t,t),mustar,sigstar)) ;
    jfr(t,1)           = probs(t,t)*lam*(1-logncdf(wstar(t,t),mustar,sigstar)) ;
    meanacceptedw(t,1) = integral(wint,wstar(t,t),Inf)./(1-logncdf(wstar(t,t),mustar,sigstar));
    expectedoff(t,1)   = lognstat(mu(t),sigstar) ;
    if t==1
        cumulsum(t,1)   = cumsum_dat(1);
        cumcomp(t,1)    = meanacceptedw(t,1) ;
    else
        cumulsum(t,1)   = cumulsum(t-1,1)+(1-cumulsum(t-1,1)).*jfr(t,1) ;
        cumcomp(t,1)    = (cumulsum(t-1,1).*cumcomp(t-1,1)+(1-cumulsum(t-1,1)).*jfr(t,1).*meanacceptedw(t,1))./((1-cumulsum(t-1,1)).*jfr(t,1)+cumulsum(t-1,1)) ;
    end
    proboffer(t,1)     = probs2(t)*lam ;
    jfrh               = jfr(t,1) ;
    meanwh             = jfr(t,1).*meanacceptedw(t,1) ;                  
    for h=t+1:capT2
        jfrh    = jfrh+(1-probs(h-1,t)*lam*(1-logncdf(wstar(h-1,t),mustar,sigstar))).*probs(h,t).*lam*(1-logncdf(wstar(h,t),mustar,sigstar)) ;
        meanwh  = meanwh+(1-probs(h-1,t)*lam*(1-logncdf(wstar(h-1,t),mustar,sigstar))).*probs(h,t).*lam*(1-logncdf(wstar(h,t),mustar,sigstar)).*integral(wint,wstar(h,t),Inf)./(1-logncdf(wstar(h,t),mustar,sigstar)) ;
    end
    expectedw(t,1)     = meanwh./jfrh ;
    %
end
fdist = jfr.*cumprod((1-jfr))./(1-jfr) ;
fdist = fdist./sum(fdist) ;
rej_atleast1 = sum(fdist.*(1-cumprod((1-rej))));

temp1 = fitlm(linspace(1,capT+1,capT+1),cumcomp(1:capT+1)) ;
temp1.Coefficients.Estimate(2);

temp2 = fitlm(linspace(1,capT+1,capT+1),cumcomp_dat(1:capT+1));
temp2.Coefficients.Estimate(2) ;


fval = 100*[(cumcomp(1)-cumcomp_dat(1))/cumcomp_dat(1) (temp1.Coefficients.Estimate(2)-temp2.Coefficients.Estimate(2))/temp2.Coefficients.Estimate(2) ...
            (expectedw(2)-beliefs_dat(1))./beliefs_dat(1) (expectedw(8)-beliefs_dat(2))./beliefs_dat(2) ...
                (cumulsum(capT+1)-cumsum_dat(capT+1))/cumsum_dat(capT+1)] ;

if isnan(fval)
    display('stop')
end